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INTRODUCTION 

Unguided  fin-stabilized  vehicles  find  wide  applications 
in  aircraft  ordnance  and  sounding  rockets  due  to  their 
relative  simplicity  and  low  cost. 

One  of  the  more  important  problems  encountered  in  un¬ 
guided  vehicle  design  is  that  both  the  nutational  and  roll 
frequencies  may  be  equal  at  some  time  during  the  flight. 

For  certain  combinations  of  vehicle  and  flight  parameters 
the  ensuing  pitch^-roll  resonance  condition  can  cause 
sufficient  amplification  of  the  vehicle  angle  of  attack  to 
produce  structural  failure  or  significant  deviation  from  the 
intended  flight  path.  The  designer  must  attempt  to  ensure 
that  pitch- roll  resonance  is  not  encountered.  In  the  event 
that  this  is  not.  possible,  he  must  ensure  that  the  crossover 
of  the  two  frequencies  takes  place  as  rapidly  as  possible • 

When  the  dynamic  pressure  and  spin  rate  are  nearly 
constan  ,  and  the  aerodynamic  forces  and  moments  are  linear 
(i.e.  no  dependence  on  roll  angle  and  linear  dependence  on 
angle  of  attack)  then  analytic  approaches  to  this  problem 
are  quite  satisfactory  (ref.  (1)),  If  these  rather  strin¬ 
gent  conditions  are  not  fulfilled,  then  a  six-degree-of- 
freedom  method  of  analysis  will  probably  be  required  for  an 
accurate  assessment  of  the  problem.  The  primary  purpose  of 
this  report  is  to  present  a  computer  program  which  can  be 
used  for  this  type  of  analysis  as  well  as  the  analysis  of 
spin-stabilized  vehicles. 

This  program  is  similar  to  those  described  in  refer¬ 
ences  (2)  and  (3)  because,  among  other  things,  the  Euler 
angle  formulation  of  the  problem  is  not  used.  The  reason 
for  this  is  twofold.  First,  each  angle  occurring  in  the 
initial  construction  of  the  inert ial-to-body-axes-transfer 
matrix  can  be  readily  defined  in  terms  of  the  physical 
variables  describing  the  motion  of  the  vehicle.  Second,  the 
singularity  occurring  when  one  of  Euler’s  angles  is  equal  to 
90°  is  eliminated. 

There  are  several  special  features  which  distinguish  this 
program  from  many  other  six-degree-of-freedom  trajectory 
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prograss  written  in  the  past. 

a.  Aerodynamic  force  and  moment  coefficients  which  are 
arbitrary  functions  of  Mach  number,  total  angle  of  attack  and 
roll  angle  are  introduced  in  table  fora.  Correct  use  of  the 
options  afforded  in  the  handling  of  the  tabulated  data  can 
result  in  considerable  time  savings  in  preparation  of  the  tables 
and  running  time  of  the  program. 

b.  The  integration  subroutine  provides  the  user  with  the 
option  of  automatically  adjusting  the  integration  time  step  on 
the  basis  of  estimates  of  the  truncation  error.  This  feature 
is  quite  important  when  the  dynamic  pressure  or  spin  rate  vary 
over  wide  ranges  during  the  vehicle  flight. 

c.  Nonrolling  body  axes  can  be  used  even  when  aerodynamic 
asymmetries  exist.  This  is  important  when  the  vehicle  is 
rapidly  spinning  since  considerable  time  savings  over  using 
rolling  body  axes  can  be  realized. 

Along  with  the  mathematical  description  of  the  program,  a 
program  listing  along  with  Identification  of  the  key  arrays  is 
included.  The  reasons  for  doing  this  are  to  insure  that  the 
user  has  complete  understanding  of  the  program  and  facilitate 
any  changes  which  may  be  desired. 

SYMBOLS 

body  reference  area 

azimuth  of  in  local  frame  (Figure  3) 

A 

azimuth  of  ?E  in  local  frame  (Figure  3) 
azimuth  of  local  wind  (Figure  3) 
speed  of  sound 

aerodynamic  force  coefficients  defined  xn  rolling  or 
nonrolling  body  frame  (Equation  4) 

aerodynamic  moment  coefficients  defined  in  rolling 
or  nonrolling  body  frame  (Equation  5) 

aerodynamic  force  coefficients  defined  in  aero¬ 
dynamic  data  frame 

aerodynamic  moment  coefficients  defined  in  aero¬ 
dynamic  data  frame 


a. 

E 

Nr 

c 


CX,CY’CZ 

CL’CM’CN 


c  ,c  ,c 

x’  y*  z 


C,,C  ,C 
V  m*  n 
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dC 


x 


body  reference  diameter 

atmospheric  density  deviation  function 
(Equation  21) 

aerodynamic  tf  rces  defined  in  rolling  or  nonrolling 
body  frame  (Equation  4) 

acceleration  due  to  gravity  at  earth's  surface 

altitude  of  body  above  earth’s  surface 

axial  moment  of  inertia  of  vehicle 

transverse  moment  of  inertia  of  vehicle  about 
body  y  axis 

transverse  moment  of  inertia  of  vehicle  about 
body  z  axis 
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I  ■ l  inertial  to  earth  frame  transfer  matrix 

(Equation  33) 

body  to  inertial  frame  transfer  matrix 
(Equations  1  and  28) 


inertial  to  local  frame  transfer  matrix 
(Equation  13) 

elements  of 


m 


mass  of  vehicle 


M  Mach  number  (Equation  19) 

M  ,M  ,M  aerodynamic  moments  defined  in  rolling  or  non- 

x  y  z  rolling  body  frame  (Equation  5) 


Q 


components  of  angular  velocity  of  vehicle  in 
body  frame 

dynamic  pressure 


radius  of  earth 

longitudinal  range  with  respect  to  the  earth 
(Equation  29) 


latitudinal  range  with  respect  to  the  earth 
(Equation  30) 


t  time 

VA  velocity  of  vehicle  with  respect  to  local  air 

mass 


speed  of  vehicle  with  respect  to  local  air 

nass  "  Pa  I 


velocity  of  vehicle  with  respect  to  earth 
speed  of  vehicle  with  respect  to  earth  - 

speed  of  local  wind 
components  of  V.  in  body  frame 


V  ,V  ,V  components  of 
XLA  yLA  ZLA 


in  local  frame 
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V  ,V  ,V 
XLE  *LE  ZLE 


W 

W 


w 


A 

B 

E 


>vl,yl,z 

w 


L 

R 


components  of  V  in  local  frame 

£ 

aerodynamic  data  frame  axes 

rolling  or  nonrolling  body  frame  axes 

earth  frame  axes 

local  frame  axes 

inertial  frame  axes 


angle  of  attack  (Equation  7) 
total  angle  of  attack  (Equation  9) 
angle  of  sideslip  (Equation  8) 

elevation  angle  of  in  local  frame  (Equation  27) 
elevation  angle  of  V_  in  local  frame  (Equation  24) 

h 

average  fin  cant  angle 

effective  configurational  asymmetry  angle  in 
pitch  plane 

effective  configurational  asymmetry  angle  in 
yaw  plane 

orientation  of  plane  of  total  angle  of  attack 
(a)  and  plane  of  elevation  of  (y^)  (Figure  4) 


P 


atmospheric  density 


T 


R 


o 


t 


0 


standard  atmospheric  density 

earth  longitude  (Equation  31) 

inertial  longitude  (Equation  14) 

orientation  of  rolling  body  frame  with  respect 
to  X£-ZA  Plane  (Equation  12) 

orientation  of  rolling  body  frame  with  respect 
to  nonrolling  body  frame  (Equation  11) 

orientation  of  rolling  or  nonrolling  body 
frame  with  respect  to  XA"ZA  plane  (Equation  10) 
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aerodynamic  roll  angle 
earth  latitude  (Equation  32) 
inertial  latitude  (Equation  15) 
angular  velocity  of  earth 

iL 

dt 

ill 

dt* 

transpose  of  matrix 
inverse  of  matrix 


(  )  value  of  variable  at  which  calculations  are 

automatically  terminated 


MATHEMATICAL  FORMULATION  OF  PROBLEM 


The  equations  of  motion  are  written  in  a  form  such  that  the 
force  equations  are  integrated  in  the  inertial  frame  and  the 
moment  equations  are  integrated  in  the  body  frame.  They  are 
also  written  in  a  form  which  allows  either  a  rolling  or  non¬ 
rolling  body  frame  to  be  used. 


* 


» 


The  force  equations  in  the  inertial  frame  are 
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and  the  moment  equations  in  the  inertial  frame  are 


Ip 

I  P 

M 

X 

X 

X 

I  q 

y. 

--  [°] 

i  q 
y 

+ 

M 

y 

I  r 

I  r 

M 

z 

z 

z 

where 


q  kp  0 


(2) 


If  the  rolling  body  axes  are  used  then  k«l  and  if  the  non¬ 
rolling  body  axes  are  used  k-0.  When  I  /  I  ,  then  the 
rolling  body  axes  must  be  used.  y 


The  derivatives  of  the  elements  of  the  direction  cosine 
matrix  [^B  ]  are  given  by  the  expressions 


ij 


3 

£ 

k-1 


ik 


% 


ij  -  1,2,3 


(3) 


Equations  1,  2  and  3  form  a  set  of  18  first  order  simul¬ 
taneous  differential  equations  describing  the  motion  of  the 
vehicle. 

Aerodynamic  Force  and  Moment  Equations 


The  aerodynamic  forces  and  moments  in  the  body  frame  are 
now  defined.  The  definitions  apply  for  both  rolling  and  non¬ 
rolling  body  frames.  The  forces  are 


(4) 
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and  the  moments  are 


M 

Cr 

X 

L 

M 

y 

-  §Ad 

c« 

M 

C 

z 

N 

(5) 


The  force  and  moment  coefficients  appearing  in  the  above 
equations  are  now  written  in  terms  of  force  and  moment 
coefficients  defined  in  the  aerodynamic  data  frame ,  with  the 
exception  of  the  moments  due  to  configurational  asymmetries 
which  are  defined  in  the  rolling  body  frame.  The  aerodynamic 
data  and  body  axes  systems  are  shown  in  Figure  1.  The 
expressions  for  the  coefficients  appearing  in  equations  (4)  and 
(5)  are 


Cv  -  C  +  C 

XXX  . 

p  A 


C_  -  C  cos  0  +  C  sin  0 
Y  n  z 


C_  -  C  cos  0  -  C  sin  0 
Z  z  y 

C_  «  C,  +  C  5  +  C, 

L  t  v  2V. 

p  A 


C.  -  C  cos  0  +  C  sin  0  +  C 
Mm  n 


Pd 


sin  0 


n  2V. 

d  P  A 

+  C  ts?  +  C  ($.  cos  <p'  -  5  sin  <p#) 
m  Av.  m*  a  d 

q  A  6a 

C„  -  C  cos  0  -  C  sin  0  +  C  •Jtk  cos  0 
N  n  m  n  2V. 

P  A 

+  C  +  C  (5_  cos  (p 9  +  6.  sin  w') 

m  2V.  m  '  B  A 

q  A  6a 


(6) 


where  C.C  .  C.  ,  C,  ,  C  ,  C  and  C  are  functions  of 
x  x  -c  '  i,  m  n  n. 

P  5  P  Q  P 

M  and  a  and  C  ,  C  ,  C.}  C  and  C  are  functions  of  M,  5  and  0'. 
y'  zf  V  m  n 


8 
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When  looking  up  C  appearing  in  the  expression  for  Cu,  ct  is' 

in  M 

q 

taken  to  b©  |  a  |  and  in  the  equation  for  C„,  a  is  taken  to  be 

1*1  •  N 


Aerodynamic  Orientation  Angles 


In  this  section  the  aerodynamic  orientation  angles,  re¬ 
quired  for  the  calculation  of  the  force  and  moment  coefficients 
given  in  equation  (6),  are  defined.  The  angle  of  attack, 
sideslip  and  total  angle  of  attack  are  shown  in  Figure  1.  It 
should  be  noted  that  for  a  given  orientation  of  the  body  with 
respect  to  the  velocity  vector  the  values  of  a  and  ft  depend 
upon  whether  the  body  axes  are  rolling  er  nonrolling.  The 


value  of  a,  however,  is  independent  of  the  choice  of  body  axes. 
Also  note  that  VA  and  a  always  lie  in  the  plane. 


a  -  tan 


+180®  *  a  >  -180° 


V 


+180°  *  0  >  -180° 


(8) 


NOLTR  64-225 


a  -  tan 


-1 


a 

BA  + 


180°  *  5  *  0° 


(9) 


i  - 

Range  of  a 

T 

~ — 

*BA 

r-v - 

yBA 

y 

ZBA 

90°  >  5  >  0° 

+ 

180°  >  5  >  0° 

] 

i 

_ 

5-0° 

i 

0 

0 

0 

The  roll  angle  0,  as  is  the  case  with  a  and  ft,  depends 
upon  whether  or  not  the  body  frame  is  rolling 


0  - 


360°  >0*0° 


(10) 


Range  of  0 

V 

yBA 

V 

ZBA 

90°  >  0  >  0° 

*  1  i 

+ 

+ 

180°  >  0  >  90° 

+ 

— 

270°  >  0  >  180° 

360°  >  0  >  270° 

o 

O 

1 

•<& 

0 

0 

Initially  the  rolling  and  nonrolling  body  frames  are 
coincident.  If  it  is  desired  to  use  a  nonrolling  body  frame 
when  the  vehicle  has  a  configurational  asymmetry  or  other  roll 
dependent  forces  and  moments,  then  the  roll  orientation  of  the 
rolling  body  frame  with  respect  to  nonrolling  body  frame  Op') 
must  be  calculated.  This  roll  orientation  angle  is  given  by 
the  expression 

<p'  -  (l-k)p  360°  >  <p'  *  0°  (11) 
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where  again 

^  ^  0  nonrolling  body  axes 

1  rolling  body  axes 

The  roll  orientation  of  the  rolling  body  axes  with  respect 
to  the  plane  of  the  total  angle  of  attack  (p)  is  given  by 

o  -  0  +  o'  (12) 


In  the  case  of  cruciform  finned  bodies  the  roll  dependent 
aerodynamic  coefficients  can  have  a  90°  or  180°  symmetry  in  <p. 
If  such  a  symmetry  exists,  the  amount  of  aerodynamic  data 
required  can  be  reduced  considerably  by  making  the  coefficients 
functions  of  the  dummy  variable  0'  instead  of  <p. 

If  there  is  a  90°  symmetry  in  the  roll  dependent 
coefficients  then 

0'  -  <p  90°  *  <p  *  0° 

0'  »  to  -  90°  180*  *  <p  >  90° 

0'  m  <p  -  180°  270°  *  tp  >  180° 

0'  -  <p  -  270°  360°  >  <p  >  270° 

and  if  there  is  a  180°  symmetry  in  the  roll  dependent 
coefficients  then 

0'  -  (p  180°  *  <p  *  0° 

0'  -  <p  -  180°  360°  >  (0  >  180° 

For  all  other  cases  0'  -  <p  and  the  roll  dependent  data 
must  be  supplied  over  the  range  360°  >  0#  *  0°. 

Functional  Calculations 

The  functional  calculations  are  those  which  roust  be 
performed  in  order  for  the  program  to  operate.  The  aerodynamic 
force  and  moment  calculations  described  in  the  previous  section 
are  of  course  program  calculations  but  are  described  separately 
because  of  their  significance.  The  remainder  of  the  functional 
calculations  described  in  this  section  essentially  express  the 
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variables  upon  which  the  aerodynamic  forces  and  moments  depend, 
in  terms  of  the  integrated  variables.  Specifically  the 
variables  Q,  M#  VA,  V  ,  V  ,  and  V  must  be  defined  in 

A  *BA  yBA  ZBA 

terns  of  XR,  YR,  ZR,  iR,  fR,  ZR,  and  . 

The  components  of  the  velocity  of  the  body,  with  respect 
to  the  earth,  in  the  local  axes  (L)  are  given  by  the  expression 


(13) 


The  direction  cosine  matrix 
expression 


j" lr  is  given  by  t;  o 
L 


R 


c 


R 


0 


-s .  s 
ip  7 
R  R 

c ,  s 

Tn 


R 


The  angles  T_  and  ipR  are  shown  in  Figure 
as  follows 


2  and  are  defined 


«  tan 


360°  >  T  *  0°  (14) 

K 


Range  of  rR 

2r 

yr 

90°  >  T_  >  0° 

R 

+ 

+ 

180°  >  T  >  90° 

K 

+ 

- 

270°  >  T  >  180° 

K 

- 

- 

360°  >  T_  >  270° 

K 

- 

+ 

tr  •  0 

6 

0 
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+180d*  tf>R  >  -180°  (15) 


Range  of  $>R 

XR 

180°  >  0  >  0° 

R 

+ 

o °  >  K  >  -i8^0 

K 

*R  “  ° 

0 

The  components  of  the  velocity  of  the  vehicle,  with 
respect  to  the  local  air  mass,  in  the  local  frame  are  given 
by  the  expression 


V 

v 

vw 

cos  Aw 

LA 

LE 

V 

Y 

yLA 

- 

yLE 

V 

nr 

V 

7  j. 

vw 

sin  S 

LAJ 

LE 

(16) 


The  winds  are  assumed  to  be  locally  horizontal  with  no  vertical 
component  (see  Figure  3) .  In  the  absence  of  winds  the 
atmosphere  is  assumed  to  rotate  with  the  earth  and  be  uniform 
in  density  for  a  given  altitude.  The  local  wind  speed  (V  ) 
and  azimuth  (Atf)  are  measured  with  respect  to  the  earth.  WThe 
wind  azimuth  has  the  range  360°  >  i  0°  and  is  the  direction 
from  north  from  which  the  wind  is  moving. 

The  components  of  the  velocity  of  the  body  with  respect 
to  the  local  air  mass  in  the  body  (B)  frame  are  given  by 


V 

V 

x_ . 

BA 

XLA 

V 

r.  -,T  r  ,  -i  T 

V 

yBA 

-  [  b]  [  'l] 

yLA 

V 

V 

_ZBA_ 

ZLA 

13 


(17) 


NOLTR  64-225 


The  elements  of  the  direction  cosine  matrix  J  are  the  terms 
^ij  ol3taine<i  the  integration  of  the  equations  of  motion. 

The  total  velocity  of  the  body  with  respect  to  the  local 
air  mass  is 


V.  «  \/v  3  +  V  a  +  V  * 

A  V  XBA  yBA  ZBA 


(18) 


The  Mach  number  is  given  by  the  expression 


M-I* 

c 

where  the  speed  of  sound  (c)  is  given  in  reference  (4) . 


(19) 


The  altitude  of  the  vehicle  is  given  by  the  expression 


h  -  yv +  v +  v  -  re 

and  the  dynamic  pressure  is 


(20) 


Q  -  iPVAa 

The  atmospheric  density  is  calculated  by  the  expression 

p  -  PA( D  +  1)  (21) 

where  p.  is  the  standard  atmospheric  density  given  in  refer¬ 
ence  (4J  and  D  is  the  density  deviation,  given  as  a  function 
of  altitude,  which  allows  the  use  of  an  arbitrary  density 
profile. 

Initial  Calculations 


These  calculations  must  be  performed  in  order  to  convert 
the  input  data  into  the  variables  required  for  the  initiation 
of  the  integration  process.  This  means  that  XR,  YR,  ZR,  XR,  fR, 

Z  ,  and  l  must  be  expressed  in  terms  of  h,  V  A  ,  y„f  r_, 
R  lj  E  E  E  E  E 
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0,  a,  0,  p,  q,  and  r. 

Since  the  earth  and  interial  frames  are  initially  coin¬ 
cident,  the  initial  values  of  the  inertial  longitude  and 
latitude  are 


The  ingle  of  attack  and  sideslip  angle  are 

a  -  a  cos  0 
fl  -  5  sin  0 

The  initial  values  of  the  inertial  coordinates  of  the  body 
(XR,  Y  ,  and  Z  )  are  calculated  using  equations  (14),  (15)  and 
the  equation 

h  "  A*  +  V  +  ZR2  -  RE  (22> 

The  components  of  VE  in  the  local  frame  can  now  be  calcu¬ 
lated  from  the  expressions  for  the  earth  referenced  flight 
path  azimuth  and  elevation  angles 


v  ! 

A  -1  ZLE 

E  “  ta”  V 

L  XLE  J 

360°  >  A£ 

Range  of  AE 

V 

ZLE 

V 

XLE 

90°  >  A_  >  0° 

E 

+ 

+ 

+ 

+ 

270°  >  A  >  180° 

£ 

- 

+ 

360°  >  A  >  270° 

£ 

- 

+ 

AE  -  0 

0 

'  0 

s  0° 


(23) 
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-90°  *  y£  *  +90° 


Range  of  y_ 

Jbi 

V 

yLE 

90°  *  yE  <  o 

+ 

0 0  >  V  *  -90° 

IS 

- 

o 

i 

w 

0 

(24) 


along  with  the  equation 


(25) 


Having  the  components  of  V£  in  the  local  frame,  the  initial 
values  of  the  inertial  velocity  components  in  the  inertial 
frame  can  be  obtained  by  inverting  equation  (13) . 

The  initial  elements  of  the  direction  cosine  matrix 
must  now  be  calculated.  The  matrix  J  is  the  trans¬ 
pose  of  the  direction  cosine  matrix  formed  by  rotating  an  axes 
system,  initially  coincident  with  the  inertial  axes,  through 
a  series  of  angles  so  that  in  its  final  orientation  it  is 
coincident  with  the  body  frame.  The  sequence  of  rotation  is 

T  about  X 
K 

-$  about  Z 
R 

-A.  about  Y 
A 

y .  about  Z 
A 

90*-©  about  X 


we 
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a  about  T 
0  about  X 

The  inertial  longitude  rR  and  latitude  $R  shown  in  Figure  2 
are  giver  by  equations  (14)  and  (15). 

The  azimuth  and  elevation  angles  of  the  vehicle's  velocity 
vector  with  respect  to  the  local  air  mass  are  shown  in 
Figure  3  and  are  defined  by  the  equations 


360°  >  A.  *  0°  (26) 

A 


90°  *  y.  *  +90°  (27) 

A 


The  »xgn  conventions  for  and  y  ^  are  the  same  as  those  for 

Ar  and  yE  given  in  equations  (23)  and  (24) .  The  components 
VE  in  the  local  frame  are  calculated  from  equations  (23),  (24), 
and  (25) .  These  components  used  in  equation  (16)  give  the 
velocity  components  V  ,  V  ,  and  V  required  in  equations 

LA  yLA_  ZLA 

(26)  and  (27) .  The  angles  6,  5,  and  0  are  shown  in  Figures  1 
and  4,  and  are  input  data.  The  angle  ©  is  measured  in  a 
plane  perpendicular  to  V.  and  for  ©  -  0°  or  180°  XA  lies  in 
the  plane  of  v  The  initial  elements  of  \  (  1  are  obtained 
by  expanding  tne  following  matrix  L 
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COS  0R 

-sin  j|)r 

0 

1 

0 

0 

sin  j|jr 

cos  0R 

0 

0 

cos  rR 

sin  tr 

0 

0 

1 

0 

-sin  rR 

COS  T 

K 

— 

Auxiliary  Calculations 


(28) 


Calculations  which  are  optional  are  described  in  this 
section. 

The  longitudinal  and  latudinal  ranges  are  given  by  the 
expressions 


R 


R_T 
£  E 


(29) 


and 


where 


R 


¥ 


“  Ve 


E 


tan 


-1 


1 

r 

E 


ipE  -  tan 


-1 


E 


\/VTV 


360°  5s  T  0l> 
E 


+180°  5>  </>_  >  -180° 
E 


(30) 


(31) 


(32) 


The  sign  conventions  on  the  earth  longitude  and  latitude  are 
the  same  as  those  in  equations  (14)  and  (15).  The  longi¬ 
tudinal  range  is  always  positive  and  the  latitudinal  range 
has  the  sign  of 

The  earth  position  coordinates  of  the  body  are  obtained 
from  the  inertial  coordinates  by  the  expression 


x„ 

E 

»—  — I 

R 

Y„ 

«  k 

Y_. 

E 

L  Ej 

R 

z„ 

Zn 

E 

R 

(33) 
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noting  that  initially  the  earth  and  inertial  franes  are 
coincident  (see  Figure  5) .  The  direction  cosine  matrix  F l 
is  given  by 


0 

sin  d^t 
cos  u^t 

The  time  rate  of  change  of  the  aerodynamic  roll  angle  is 


a 


r  f 

Ff 

V 

y  v 

V 

Z  V 

ZBA 

— * - r  x_. 

L  m  BA_ 

yBA 

L~ +  q  v 

V  8  V  * 


+ 


+  p 


(34) 


DESCRIPTION  OF  MACHINE  PROGRAM 

This  program  is  coded  in  FORTRAN  IV  and  has  been  prepared 
for  running  on  an  IBM  7090  with  the  IBSYS  monitor.  It 
consists  of  a  routine,  entitled  MAIN,  which  exercises  overall 
control  of  the  program  and  fourteen  subroutines.  Descriptions 
of  all  subroutines  are  given  in  Appendix  A  and  FORTRAN  listings 
of  MAIN  and  the  subroutines  SETUP,  AERO,  DERIV,  TERM,  OUT  and 
FNOL-3  are  given  in  Tables  1  through  7. 

The  most  important  arrays  appearing  in  this  program  are 
the  KY,  Y,  C,  and  D  arrays.  These  arrays  are  identified  in 
Tables  8  through  11.  The  KY’s  are  program  options  and  controls 
and  identify  tabulated  data  described  below.  The  Y  array 
contains  additional  program  controls  and  most  of  the  program 
variables  appearing  in  the  equations  of  notion  and  the  D's 
are  their  first  derivatives. 

Input 


The  cards  necessary  to  run  this  program  contain  controls, 
initial  conditions,  constants  and  tabulated  data. 

The  tabulated  data  consist  of  aerodynamic  coefficients, 
the  density  deviation  (D)  and  the  wind  velocity  (V^)  and 
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wind  azimuth  (A^) .  The  aerodynamic  coefficients  depend  upon 
Mach  number,  angle  of  attack  and  roll  angle  as  explained  in 
the  section  "Aerodynamic  Force  and  Moment  Equations."  The 
density  deviation,  wind  velocity  and  wind  azimuth  are  functions 
of  altitude  alone.  These  tables  are  prepared  in  accordance  with 
the  instructions  in  Appendix  A,  section  7.  These  data  may  be 
introduced  into  machine  memory  directly  from  cards  or  from 
cards  to  special  tape  and  then  to  machine  memory. 

a.  Tables  introduced  directly  into  memory  from  cards  - 
When  using  this  mode  of  operation  subroutine  MEMTAB  places 
the  tables  in  memory.  The  card  deck  required  in  this  case  is 
as  follows: 

Card  No.  Item  Format 

1  $JOB 

2  $ EXECUTE  IBJOB 

3  $IBJOB  NOSOURCE 

Trajectory  Program  Binary  Deck 

4  $DATA 

5  ITAB,  JO 

6,7  KT(21)  -  KY(JO) 

Tabulated  Data 

8  Blank 

9  I  NT 

10,11  KY(1)  -  KY(20) 

12-19  Y (200)  -  Y(INT) 

Repeat  Cards  9-19 
for  Each  Additional  Run  Required 

20  Blank 

21  Termination  (7-8  punch 

1st  col.) 

For  this  mode  of  running,  ITAB-0.  The  number  of  KY 
inputs  is  JO.  If  there  are  no  winds  then  the  tables  and 
A  are  not  required  and  JO-33  and  if  winds  are  present  JO-35. 

The  number  of  Y  inputs  is  INT-200.  For  this  program  INT=*238. 


1017 


1017 

If 

5E14.6 
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When  the  tables  are  introduced  into  memory  in  this 

manner  they  must  appear  in  the  following  order:  C  ,  C  ,  C  , 

C  ,  C,,  C,  ,  C,  ,  C  ,  C  ,  C  ,  C  ,  C  ,  D,  V_.,  X  xp  y 

zf  V  l'  l  1  to*  nf  m  *  n  *  m.  *  *  W*  r 

6  p  q  p  6a 

A^.  If  multiple  runs  are  made  then  all  runs  must  utilize  the 
same  tabulated  data. 

b.  Tables  placed  on  special  tape  and  then  introduced 
into  memory  -  The  tables  must  be  placed  on  tape  by  a 
separate  program,  TAPER.  When  this  is  done  the  tables  are 
then  transferred  to  machine  memory  by  subroutine  TABLIN. 

The  required  card  deck  for  this  mode  of  operation  is  as 
follows: 


Card  No. 


_ Item 

$JOB 

$EXECUTE 

$IBJOB 


Format 


TAPER 


$DATA 


Tabulated  Data 
Blank 

Termination 
$JOB 
$ EXECUTE 
$IBJOB 


Trajectory  Binary  Deck 


10 

$DATA 

11 

ITAB,  JO 

1017 

12 

INT 

tf 

13-16 

KY (1)  -  KY(JO) 

tl 

17-^4 

Y (200)  -  Y (INT) 

5E14.6 

Repeat  Cards  12-24 
for  Each  Additional  Run  Required 

25  Blank 

26  Termination  (7-8  punch 

1st  col.) 
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For  this  mo-Ja  of  running  ITAB-1.  JO  and  INT  are  as 
previously  described.  If  the  tabulated  data  are  to  be  supplied 
by  a  previously  prepared  tape  then  the  first  card  of  the  deck 
will  be  card  number  7. 

When  tables  are  placed  on  the  special  tape  they  are 
implicitly  numbered  in  the  order  in  which  they  appear.  The 
numbers  assigned  to  KY(21)  through  KY(35)  on  the  input  cards 
assign  a  given  table  on  the  special  tape  to  a  particular  KY. 

In  other  words,  if  KY(21)-1  then  the  first  table  on  the  tape 
contains  and  so  on.  It  should  also  be  noted  that  more 
tables  than  are  required  for  a  given  run  may  be  placed  on  the 
special  tape.  The  advantage  of  doing  this  can  be  seen  for  the 
case  when  a  series  of  runs  are  desired  each  having  a  different 
value  for  C  .  The  tables  would  be  prepared  just  as  they 
would  for  a  single  run  but  with  additional  tables  of  C 
included.  The  value  assigned  to  KY(21)  would  change,  from 
run  to  run,  depending  on  the  particular  table  of  C  desired. 

The  numbers  assigned  to  KY(22)  through  KY(35)  woulS  remain  the 
same.  This  feature  is  not  available  when  the  tables  are  read 
directly  into  memory  by  MEIiTAB  since,  in  this  case,  C  is 
always  the  first  table,  C  the  second,  and  so  on.  x 

p 


Options  and  Controls 

Options  and  controls,  other  than  those  which  must  be 
exercised  in  order  to  use  the  required  subroutines,  will  be 
discussed  in  this  section. 

a.  Rolling  or  nonrolling  body  axes  -  The  value  of  KY(12) 
determines  whether  the  body  axes  roll  with  the  body  or  are 
nonrolling. 

KY(12)  -  0,  nonrolling 
KY{12)  -  1,  rolling 

The  only  restriction  on  the  selection  of  the  type  of  axes  to 
be  used  is  that  if  I  r*  I  the  rolling  frame  must  be  used.  The 
aerodynamic  forces  aXd  moments  are  calculated  in  a  manner  such 
that  the  data  required  are  the  same  for  both  types  of  axes. 

The  nonrolling  axes  can  be  used  even  when  the  vehicle  has  an 
aerodynamic  asymmetry. 
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Several  factors  which  may  influence  the  choice  of  the  body 
axes  to  be  used  are  the  output  desired  and  the  running  time. 
Since  the  quantities  a,  f$,  0,  q,  r,  F  ,  F  ,  M  and  M  are 
defined  in  the  body  axes,  their  particular  values  wifi  depend 
upon  the  axes  chosen.  If  the  body  is  rapidly  spinning  the 
use  of  nonrolling  axes  will  in  general  allow  a  significantly 
larger  integration  tine  step  than  will  the  rolling  axes. 

b.  Aerodynamic  symmetry  -  This  option  applies  only  to  the 
aerodynamic  coefficients  which  are  functions  of  the  aerodynamic 
roll  angle  C  ,  C  ,  Cf,  and  Cq.  The  value  assigned  to  KY(ll) 
indicates  thC  degree  of  roll  symmetry  which  the  coefficients 
have. 


KY(ll)  -  1,  90*  symnetry 
KY(ll)  -  2,  180* 

KY(ll)  -  3,  360* 

For  the  case  of  90°  or  180°  symmetry  in  roll  the  roll  dependent 
coefficients  need  only  be  tabulated  for  roll  angles 
O°^0#i9O°  or  O®*0'*18O#,  respectively,  instead  of  the  full 
range  0°  to  360°.  This  option  can  result  in  considerable  time 
savings  in  data  preparation  and  reduction  of  storage  required 
for  aerodynamic  data.  A  more  complete  discussion  of  this  option 
is  given  in  the  section  "Aerodynamic  Orientation  Angles." 

c.  Force  free  calculations  -  Above  the  altitude  Afl 
(Y(203))  the  aerodynamic  forces  and  moments  are  set  equal  to 
zero.  This  option  considerably  shortens  the  computation 
procedure  when  this  type  of  calculation  is  desired. 

d.  Winds  -  If  there  are  no  winds  KY(7)-rq  and  if  winds 
are  present  KY(7)-1.  In  the  event  that  there  are  no  winds 
the  wind  velocity  and  azimuth  tables,  designated  by  KY(34)  and 
KY(35)  respectively,  do  not  have  to  be  included. 

e.  Rotating  earth  -  If  Y(204)  is  equal  to  zero  then  the 
earth  is  treated  as  nonrotating.  The  effects  of  a  rotating 
earth  are  included  by  setting  Y(204)«  Wg  rad/sec. 

f.  Correction  of  -  This  matrix  must  remain  very 

nearly  orthogonal  throughout  the  calculational  procedure  in 
order  to  accurately  predict  the  flight  path  of  the  body. 
Experience  has  shown  that  errors  in  this  matrix  can  be  held  to 
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a  minimum  if  it  is  corrected  3  times  each  time  step  by  the 
iterative  equation 


-1 


As  a  check  on  the  orthogon?lity  of  |  the  matrix 
can  be  printed  every  KY(10)  integration  time  steps 


g.  Output  printing  frequency  -  The  frequency  of  the  out¬ 
put  print-out  is  once  every  KY (9)  integration  time  steps. 

h.  Termination  -  Calculations  are  automatically 
terminated  when  any  of  the  conditions  below  are  satisfied. 


Output  -  The  output  data  format  is  shown  in  Table  13.  On 
the  page  preceding  the  first  page  of  output  data,  the  input 
data  are  printed. 


CONCLUDING  REMARKS 

As  previously  stated,  one  reason  for  the  detailed  docu¬ 
mentation  of  this  program  is  to  enable  the  user  to  modify  it 
to  fit  his  particular  needs  with  a  minimum  of  effort.  Once 
familiar  with  the  program,  it  is  easy  to  see  how  the  effects 
of  powered  flight,  thrust  and  mass  asymmetries  or  a  guidance 
system  can  be  incorporated  with  a  modest  effort. 

It  is  felt  that  this  program  is  capable  of  very  efficient 
operation.  The  use  of  nonrolling  body  axes  and  automatic  time 
step  adjustment,  when  the  situation  favors  their  application, 
are  proven  time  savers.  The  input  tables  allow  an  accurate 
representation  of  the  aerodynamic  coefficients  with  a  minimum 
of  time  required  for  their  preparation.  It  should  also  be 
noted  that  by  holding  the  auxiliary  calculations  and  the  out¬ 
put  to  a  minimum  considerable  savings  in  running  time  can  be 
realized. 
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APPENDIX  A 

DESCRIPTION  OF  PROGRAM  SUBROUTINES 

In  this  appendix  the  1  notion  of  each  subroutine,  as  it 
applies  to  this  program,  x  briefly  described.  The  FORTRAN  II 
version  of  FNOL-3  is  documented  in  reference  (5) .  The  sub¬ 
routines  TABLIN,  MEMTAB,  FROMRA,  ARDCFT,  SINCOS,  MATVEC, 

MATINV  and  ARKTAN  are  also  documented  and  available  from  NOL. 

For  convenience  the  call  line  variables  are  also  given. 
Unless  otherwise  indicated  T  is  time,  the  independent 
variable,  C  is  the  dependent  variable  and  D  is  its  derivative. 

1.  SETUP  (T,  C,  D)  -  This  subroutine  initializes 
variables  and  along  with  MAIN  performs  initial  calculations 
required  for  the  first  integration  of  the  equations  of  motion. 

2.  AERO  (T,  C,  D)  -  The  aerodynamic  force  and  moment 
coefficients  are  calculated  by  this  subroutine. 

3.  DERIV  (T,  C,  D)  -  The  derivatives  of  the  dependent 
variables  appearing  in  the  equations  of  motion  are  calculated 
here. 

4.  TERM  (T,  C,  D,  F)  -  The  purpose  of  this  subroutine  is 
twofold.  TERM  corrects  the  elements  of  |_£gj  *n  or<*er  to  help 
maintain  its  orthogonality.  This  routine  also  provides  for 
the  termination  of  the  calculations  as  specified  by  the  input 
data.  These  functions  of  TERM  are  discussed  in  another 
section.  The  variable  F  in  the  call  line  is  the  variable 
which  determines  termination. 

5.  OUT  (T,  C,  Dt  ERROR,  N,  L,  H)  -  Auxiliary  calculations 
which  are  not  required  in  the  actual  integration  of  the 
equations  of  motion  are  performed  in  this  subroutine.  OUT  also 
provides  for  the  storage  and  print- out  of  the  data  called  for 
in  the  output  format.  H  is  the  integration  time  step  at  time 
T.  The  terms  ERROR,  N  and  L  are  defined  in  section  6. 
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6.  FNOL-3  (J,  N,  G,  L,  M,  ENE,  T,  C,  D,  DERIV,  TERM, 

OUT)  -  This  subroutine  accomplishes  the  integration 
of  the  equations  of  motion.  It  is  the  FORTRAN  IV  version  of 
FNOL-2  described  in  reference  (5).  This  subroutine  will 
numerically  integrate  a  system  of  up  to  30  simultaneous,  first 
crder  differential  equations.  The  interval  of  integration  can 
be  automatically  adjusted  to  hold  the  absolute  or  relative 
truncation  error  within  specified  bounds.  The  call  line 
variables  are  defined  as  follows: 

a.  J«KY(3)  controls  the  mode  of  integration.  If 
automatic  time  step  adjustment  is  desired  J-2.  Under  this 
inode  of  integration  the  Runge-Kutta  method  is  used  for  the 
first  four  time  steps,  and  then  the  Adams-Moulton  predictor- 
corrector  method  is  used  until  termination.  The  integration 
is  completed  with  four  iterations  again  using  the  Runge-Kutta 
method. 


b.  N-KY(4)  is  the  number  of  differential  equations  to 
be  integrated.  For  the  program  described  in  this  report, 

N“19. 

c.  G“Y(207)  is  the  initial  interval  of  integration. 

d.  L-KY(5)  is  always  equal  to  zero  in  this  program. 

e.  M=*KY(16)  is  always  equal  to  1  in  this  program. 

f.  ENE=Y(206)  controls  the  interval  of  integration. 

If  ENE-0,  the  differential  equations  are  integrated  with  a 
constant  time  step  equal  to  G.  If  ENE>0  then  the  interval  of 
integration  is  adjusted  so  that  either  the  absolute  or 
relative  truncation  error  is  maintained  within  the  bounds 

*10-ENE  3  and  si<rENE.  If  the  automatic  time  step  adjustment 
feature  is  used,  the  program  may  restart  several  times  if  G 
is  such  that  the  initial  absolute  or  relative  truncation 
errors  are  not  maintained  within  the  specified  bounds. 

g.  ERROR=Y(205)  specifies  whether  or  not  the  time 
step  interval  is  adjusted  on  the  basis  of  absolute  (TE)  or 
relative  (RE)  truncation  error.  If  ERR0R=-1,  the  adjustment 
is  on  the  basis  of  absolute  truncation  error  and  if  ERROR=0 
it  is  based  on  the  relative  truncation  error.  The  formulas 


A-2 
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for  the  errors  are 


te- 


C(p)-C(c) 


R  - 

£ 


where  C  F  and  C  '  are  given  by  the  Adams-Moulton  predictor- 
corrector  formulas.  Care  must  be  exercised  in  the  selection  of 
absolute  or  relative  truncation  error  as  a  basis  for  time  step 


and  c<o)-  oao'5), 


adjustment.  For  example,  if  C<P)^  o(lcf3)  and  0(10~5) 

then  T  -  0(10”4)  and  R  =  0(10).  In  this  case  if  relative 
truncation  error  is  used  the  time  step  interval  may  be  adjusted 
to  a  value  many  times  smaller  than  that  required  for  accurate 
calculations. 


This  subroutine  has  greater  generality  than  indicated  here. 
Reference  (5)  should  be  consulted  for  further  details. 

7.  MEMTAB  (KY,  Y)  -  The  aerodynamic  coefficients, 
atmospheric  density  deviation  and  wind  data  required  by  this 
program  are  punched  on  cards  in  tabular  form.  These  tables 
are  introduced  into  the  machine  memory  by  this  subroutine.  KY 
is  the  irray  of  table  numbers  given  on  the  input  cards  and  Y 
is  the  array  in  which  the  tables  are  stored.  In  this  program 
the  tabulated  functions  are  functions  of  1,  2  or  3  variables, 
with  each  function  in  a  separate  table.  The  tabulation  of  a 
function  of  3  variables  would  be  as  follows: 


a.  Control  card 

L,  N,  M,  nx,  n2,  . . .  n^  (FORMAT  1415)  where 
L*=0  all  cases 

N-3  number  of  independent  variables 

M-l  each  table  contains  only  1  function 

n^n^n^  numbers  of  values  of  each  independent 

variable  for  which  values  of  the  function 
are  tabulated. 


b.  Listing  of  values  of  independent  variables  for 
which  function  is  tabulated  -  On  the  first  card/cards  the  n^ 


A- 3 
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values  of  the  first  independent  variable  are  listed.  On  the 
succeeding  cards  the  n 2  and  n  values  of  the  second  and  third 
independent  variables  are  listed.  The  following  restrictions 
apply.  Values  of  2  different  independent  variables  may  not 
appear  on  the  same  card.  At  least  2  values  of  each  independent 
variable  must  be  listed,  all  values  must  be  distinct  and  must 
be  in  ascending  order.  The  format  for  these  cards  is  6E12.7. 

c.  Listing  of  values  of  dependent  function  -  The 
three  independent  variables  are  designated  N  ,  N  and  N  .  The 
number  of  values  of  each  of  these  variables  is  n^,  n*  and  n  , 
respectively.  A  block  of  data  contains  the  values  or  the 
function  for  all  values  of  N3  listed  and  one  particular  value 
of  Nj  and  Ng.  The  first  block  corresponds  to  the  first  value 
of  N_  and  Ng  listed.  The  second  block  corresponds  to  the  first 
value  of  and  the  second  value  of  Ng.  These  blocks  are 
repeated  until  a  set  of  blocks  for  the  first  value  of  and 
all  values  of  N2  have  been  presented.  Sets  for  the  remaining 
values  of  follow  until  the  table  is  completed.  As  a  check 
there  are  n*  sets,  n^  x  n„  blocks  and  x  n2  x  n_  distinct 
values  of  the  function.  The  format  for5- this  tabulation  is 
also  6E12.7. 

Subroutine  LOCOCT  is  listed  as  a  separate  subroutine  but 
is  an  integral  part  of  MEMTAB. 

8.  TABLIN  (I,  J,  KY,  Y)  -  This  subroutine  serves  the 
same  purpose  as  MEMTAB  except  that  it  reads  tables  from  tape 
into  memory.  I  is  the  number  of  the  tape  unit  on  which  the 
tape  containing  the  tabulated  data  is  mounted.  The  number  of 
tables  on  the  tape  is  J,  which  is  automatically  set  by  the 
program.  KY  and  Y  are  the  same  as  their  counterparts  in 
MEMTAB. 


9.  FROMRA  (KY,  N,  1,  Up  U£ . UN,  V,  Z)  -  This  sub¬ 

routine  extracts  data  from  the  tables  placed  in  memory  by 
.MEMTAB  or  TABLIN.  Entering  FROMRA  with  the  N  independent 
variables  XL,  Ug,  ...,  UN  it  linearly  interpolates  or 
extrapolates  2N-1  times  and  supplies  the  desired  value  of 
the  function  V.  The  table  is  again  designated  by  KY.  Z  is 
an  array  of  temporary  storage  which  must  have  the  dimension 
of  ZN  +  2N  or  greater. 
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10.  ARDCFT  (H,  P,  T,  D,  C,  G)  -  The  earth’s  atmospheric 
properties  as  presented  in  reference  (4)  are  supplied  by  this 
subroutine  up  to  an  altitude  of  10®  feet.  Entering  with  the 
altitude  H,  the  ratios  to  the  corresponding  sea  level  value  of 
pressure  (P),  temperature  (T),  density  (D),  speed  of  3ound  (C) 
and  acceleration  due  to  gravity  (G)  are  given. 

11.  S1NC0S  (A,  SA,  CA,  N)  -  This  subroutine  supplies  sin 
A-SA  or  cos  A-CA  for  the  range  0°^A^360°.  If  N-0  A  is  in 
degrees  and  if  N-l  A  is  in  radians.  Subroutines  sin  and  cos 
must  be  included  as  part  of  the  system  library. 

12.  MATVEC  (A,  B,  C,  N)  -  Products  of  matrices  of  orders 
(3,3)  and  (3,1)  are  computed  by  this  subroutine. 


N-0, 

C-AB* 

N-l, 

C-A  B* 

N-2, 

C-AB 

N-3, 

T 

C-AB 

N-4, 

T 

C-AB1 

T  T 

N-5, 

C-AB 

The  A,  B  and  C  arrays  are  stored  column-wise  starting  at  the 
left.  The  symbol  *  indicates  that  B  in  these  cases  is  a  (3,1) 
array.  In  all  other  cases  A,  B,  C  are  (3,3)  arrays. 

13.  MATINV  (A,  B,  C)  -  This  subroutine  computes  the 
transpose  ^B)  and  inverse  (C)  of  (3,3)  matrices  (A).  If  the 
determinant  of  A  is  zero,  neither  B  or  C  is  calculated,  a 
comment  is  printed  and  control  is  returned  to  the  calling 
program. 


14.  ARKTAN  (A,  B,  C,  N)  -  Arctangents,  defined  by  the 
expression  c  _  tan“*  (A/B)  are  comPuted  by  this  subroutine. 

If  N-0,  C  is  in  degrees  and  if  N-l,  C  is  in  radians.  The 
range  of  the  output  is  -90°  ^ C  £+90* .  If  A-0  B<0,  0-90°  and 
if  B>0,  C-  +  90°.  Subroutine  ATAN  must  also  be  included  as 
part  of  the  system  library. 
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TABUS  1 

MAIN  ROUTINE  -  TRAJECTORY  PROGRAM  6UUR1 
COMMON  Y.KY 

DIMENSION  C  { 25 ) #0(25) *Y(8000) *KY ( 50 ) 
EXTERNAL  DERI V *TERM, OUT 
READ(5*1000) ITAB.JO 
IF  < ITAB)1,98.1 

98  READ ($* 1000) (KY( I ) *I*21*JO) 
WRITE(6.1002) (KY< I ) #1=21* JO) 

CALL  MEM TAB (KY(21)*Y<600) ) 

1  DO  3  I  J*1 » 500 

3  Y ( IJ)*0.0 

DO  4  I J*1 .25 
C( I J)*0.0 

4  D< IJ)*0.0 

DO  5  I J*1 *20 

5  KY ( I J ) *0 
READ(5*1000K'INT 
IF ( INT ) 100*100*2 

2  I F ( ITAB)7*6.7 

6  READ(5*1000) (KY( I ) .1*1.20) 
WRITE(6*1002)KY(1) *KY ( 2 ) 

WRITE (6. 1002)  (KY(  I  )  .1*3*20) 

GO  TO  8- 

7  RE AD (5*1000) < K Y ( I ) » 1*1. JO) 
WRITEi6*1002)KY(l).KY(2) 

WRI TE( 6.1002 ) (XY( I ) *1*3. JO) 

8  READ( 5*1001) (Y(I)»I*200tINT) 

WRITE (6. 1003) (Y( I ). I *200. INT) 

WR I TE ( 6  *1004 ) 

1000  FORMAT  (1017) 

1001  FORMAT  (5E14.6) 

1002  FOPMAT( 1H0.10I10) 

1003  FORMATt 1H0.1P5E17.6) 

1004  FORMAT  (1H1) 

I F ( ITAB)26, 22*26 

26  I F  <  K Y  <  7 ) J25.24.25 

24  KYZ*13 
GO  TO  21 

25  KYZ*15 

21  CALL  TABLIN(18.KYZ»KY(21) *Y(600)  ) 

22  CALL  SI  NCOS ( Y ( 2  2  3 ) » SN  *CN* 0  ) 

Y ( 299 )— (Y<222| )*CN 
Y(290)*(Y(222) )*SN 

CALL  S I  NCOS (Y(217)*SA.CA»0) 

CALL  SI  NCOS (Y(218)»SG*CG*0) 

Y( 66 ) *Y ( 216 )*CG*CA 
Y(67)*Y(216)*SG 
Y(68)*Y(216)*CG*SA 
IF ( <Y (7 ) )  27*27*28 

27  Y ( 160 ) *Y ( 66 ) 
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Y(161)*Y(67) 

Y ( 162 ) * Y ( 68 ) 

GO  TO  29 

28  CALL  SI  NCOS (Y(228)»SA*CA*0) 

Y( 160)=Y(66)+Y(227)#CA 

Y< 161 )=Y(67) 

Y( 162)-Y(68)+Y(22?)*SA 

29  CALL  SETUP { T »C *D) 

LB  MATRIX  COMPUTED  IN  SETUP 
C(4)  =  { Y(229)+Y(215) )#Y(47) 

C ( 5  >  *  ( Y(229)+Y<215> )*Y(46)*Y152) 

C(6 )  *  ( Y(229)+Y(215))#Y(46)*Y(53) 

Y  (  37 ) *0 • 

Y ( 38 ) *0 • 

Y<  39>*Y(204)*SQRT{C(5)**2+C{6)**2) 

CALL  MATVEC  < Y ( 72 ) • Y ( 37 ) » Y ( 40 ) . 1 ) 

CALL  MATVEC (Y(72)»Y(66)*Y(37) *1) 

DO  12  1*1.3 

C< I )*Y( 1+36 )+Y ( 1^39  > 

12  C< 1+6 ) *Y (I ♦223 ) 

C< 19)*0. 

T*Y ( 208 ) 

C ( 2 1 )  *0* 

C  (  22 ) *Y ( 205 ) 

CALL  FNOL3(KY<3) .KY(4).Y<207> .KY(5> .KYJ6).Y(206) .T.C.D. 
1DERIV. TERM. OUT) 

GO  TO  1 
100  STOP 
END 
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SUBROUTINE  SETUP(T#C»D) 

TRAJECTORY  PROGRAM  6UUR1 
COMMON  Y *KY 

DIMENSION  C( 25) #D ( 25 ) • Y ( 8000 ) *<Y ( 50 ) 

Y ( 190)«Y< 223) *.0174532923 

CALL  SI  NCOS (Y(190)*Y(51)#Y(50)  #1) 

Y< 190) *Y< 222) **0174532925 

CALL  SI  NCOS (Y(190)*YC57)#V(56)*1) 

Y ( 190 )*Y( 2 19) *#0174532925 

CALL  S I  NCOS (Y(190)*Y(59)»Y(58) *1) 

Y ( 158 ) *SQRT ( Y ( 160 )**2+Y < 162 ) **2 ) 

Y( 159)«SQRT(Y(160)**2+Y(161)**2+YU62)**2) 
Y(48)*Y(160)/Y(158) 

Y(49)*Y<162)/Y(158) 

Y(54)=Y(158)/Y<159) 

Y(55)*YC161)/Y(159) 

Y ( 157 ) *Y<  221 ) 

CALL  SINCOS  < Y C 157) • Y C 47 ) ,Y<46) .0) 

Y ( 156 ) *Y ( 220 ) 

CALL  SINCOS  (Y(156) * Y ( 53 ) * Y ( 52 ) *0 ) 

Y( 1 ) *1* 

Y ( 2 ) *0* 

Y ( 3 ) *0* 

Y ( 4  )  *0* 

Y ( 5  )  *Y ( 50 ) 

Y ( 6 ) *-Y 151) 

Y ( 7 ) *0* 

Y «  8 ) *Y  C  51 ) 

Y ( 9 ) *Y ( 50 ) 

Y ( 10 ) *Y ( 56  I 
Y ( 11)«0« 

Y ( 12 1 *Y ( 57 1 
Y ( 13 ) *0* 

Y ( 14 ) *1  * 

Y ( 15 ) *0* 

Y ( 16 ) *-Y ( 57 ) 

Y ( 1 7 ) *0 • 

Y ( 18 ) »Y ( 56 ) 

Y ( 191*1* 

Y ( 20 ) *0* 

Y( 2 1 ) *0# 

Y ( 22 ) *0* 

Y( 23 ) *Y ( 59 ) 

Y ( 24 ) «- Y ( 58 ) 

Y <  25 ) *0* 

Y ( 26 ) »Y ( 58 ) 

Y  (  27 ) *Y ( 59  ) 

Y ( 28 ) aY ( 54 ) 

Y (  29 ) «-Y ( 55 ) 

Y (  30 ) *0* 

Y  (  31 ) *Y ( 55  ) 
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Y  <  32 ) *Y ( 54 ) 

Y ( 33 ) *0» 

Y ( 34 ) *0 • 

Y ( 35 ) *0# 

Y ( 36  )  «1 • 

Y ( 163 ) *Y ( 48 ) 

Y ( 164 )  *0  • 

Y ( 165 ) *~Y ( 49 ) 

Y ( 166 ) *0« 

Y ( 167 ) *  1 • 

Y ( 168 ) *0* 

Y ( 169 ) *  Y ( 49 ) 

Y( 170)*0# 

Y  ( 171 ) *Y ( 48 ) 

Y ( 172 ) *Y ( 46 ) 

Y( 173)«Y<47) 

Y<174)*0. 

Y ( 175 ) *-Y ( 47 ) 

Y ( 1 76 ) *  Y ( 46 ) 

Y ( 177 ) *0* 

Y(178)*0. 

Y ( 179 ) *0« 

Y ( 180  >  *  1 • 

Y  (  1 8 1 )  *  1  • 

Y ( 182 ) *0« 

Y  (  183 ) *0 • 

YdSHlsQ. 

Y ( 165 ) c Y ( 52  ) 

Y ( loo ) s-Y (  53 ) 

Y ( 1 67 ) *0 • 

Y( 18d)=Y(5.3) 

Y ( 189 ) = Y ( 52 ) 

CALL  MATVEU Y(163)tY( 172) *Y(37),2) 
CALL  MATVtC(Y(28) *  Y ( 3  7 ) »Y(37) *2) 
CALL  MATVtC(Y(19)  *  Y  (  3  7  )  *Y(37)  *2) 
CaLL  maTVcL(Y(1u)  »Y  (  3  7  )  *  Y  (37)  *  2  ) 
Call  maTVcl(Y(1b1)*Y(j7)*Y(:>7)  *2 ) 
CALL  MaTVlC(Y( ^7 ) *r (1 ) *U lu ) *5) 

Y ( 72 ) *Y ( 40 ) 

Y ( 7  3 ) = Y ( 47 ) 

Y( 74)=u* 

Y(75)=-Y(47)*Y(52 ) 

Y(7o)  =  rUc)*Y(3d 
Y(77)S-Y|33) 

Y ( 7o ) =  -Y ( 47  )  *Y ( 01 ) 

Y ( 79 ) = Y ( 4b ) *Y ( PJ ) 

Y(6U)=Y C?2) 

RETuRi'* 

END 
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SUdROuTINE  AERu(T*C»u) 

TRAJECTORY  PkOokari  bOUKl 
CurihON  Y*kY 

ul  rici'ioI  un  C(25)*u(23),Y(bouu)  *f>Y(3<j)*2(:>0) 

I F ( T )  I00*l00#2uu 
iuu  J=0 

2u u  IF ( mdS( i ( 91 ) )-lu,t-7) 1  ♦  1 »4 

1  Y(91)=u. 

Y ( 1 29 ) =0 • 

I F ( Y(92) )  2  *  3  #  3 

2  Y(22^)=io0* 

Y ( 13o ) =-l • 

GO  Tu  3 

3  Y < 223 ) =  u, 

Y( 13u  )  =  1 *0 
GO  TO  3 

4  CALL  aRkTAnI Y ( 91 ) *  Y  <92 ) »Y(223) *0) 

3  Y( 19u)=oGRT( Y(9l)**2+Y(92)**2) 

C«LL  ARkTaM Y( 19u) »  Y ( 9u ) *  Y ( 222 ) » 0 ) 

Call  arkTahI Y ( 92 ) » i ( 90 ) » Y <  299 ) *0 ) 
laLL  arnTain  (  Y  (  91 )  »Y{90)  » Y ( 29b  )  #0 ) 

IF  (  Y(  cc-z  )  )  i  JfO»o 
13  Y  (  )  =  Y  (  22j  )  +joo. 

6  CALL  SIimCuM  Y(223  )  .  Y(  129)  *Y(lJ0)  *0) 

Y  (  151)=AriuD(C(19)  ,to.2b3 18331 ) *57. 2937795 
CALL  SISCOS  <  Y ( 1 3 1 ) • Y ( 191 ) • Y ( 1 92 ) *0 ) 
Y(1h9)=Y(223)+Y(131) 

IF(KY(ll)-2)  20  »22  t  25 
20  Y ( 131 ) =AriOl) ( Y ( 149 ) » 90,  ) 

GO  TO  30 

2  2  Y  (  131 )  =a‘i»iOu (  Y  (  1h9  )  *  loo* 0 ) 

60  TO  30 

2  5  Y ( 1 3 1 ) =AhOO ( Y ( 149 ) * 3b0,  ) 

30  CONTINUE 
v U  40  I=l»3 

4u  Y  ( I +1 32 ) =C ( I+6)*Y(235)*,3/Y<93) 

Call.  Frurira  (nY(21)*2*1*Y(11o)*Y(2c2)»Y(132)*2.) 

CALL  FRuriRa(kY( cc  )  *  2 » 1 > Y  (  1 lu >  t  Y <  222  )  * Y (  1 1  7  )  .2  ) 

CAll.  Fkurirm  (fvY  (<i3)  *^*1*Y(  jlio)  )  »Y  (x^l )  ♦Y' 1 13^)  *<2.) 

CALL  FROrirm  (  nY  (  24  )  *  3  *  1  » Y  (  1 lo ) * Y ( 222 ) * Y ( 1 3 1 ) t Y 1 1 34 ) *  2  ) 
CALL  FruriRA  (  nY  (25)  »3,1*Y(  110)  *Y  (  222  )  »  Y(  l3l )  tY  ( 133)  tL) 
CALL  F  aUriRA ( nY <2o)*2*i*Y(llu)>Yl222)»Y(iib)»2) 

CALi.  F  aurira  (kY(27).2,1,Y(  xlu)  *Y(222)  *Y(l3t>)*2) 

CALL  FauriRm (nY(2o)».?»1»Y{11o)*Y(2<ci)»Y(1.?1)»Y(137)»2.) 
Call  F  rui*iRA  (  sY  (  29  )  »  3  *  1  ,  Y  (  1  lo )  *  Y  (  222  )  t  Y  (  iil  )  '(l30)»2) 

Y ( lu9 ) =AD3 ( Y ( 299 ) ) 

LALL  t  AURIKA  (NY(3U)*2*l»Y(iiU)*Y(i09)#Y(i39)»2) 

Y ( luv ) = ado ( Y ( 290 )  ) 

CALL  F  aOrikA  ( &.Y  (  30  )  *2»ltY(llU)tY(l09)*Y(l40>*2) 

CALL  FROriRA ( kY ( 31 )*2»l*Y(llu) *Y(222)*Y(119)#2) 

Call  FRuhRA ( kY ( 32 >#2tl*Y(lio)*Y( 222  )  ♦Y(l03)  *2  ) 
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Y(lll)*Y(l*2)+Y(117)*Y(l33) 

Y(112)*Y(1.33)*Y(1.3U)*t(l.a4)*Y(1^9) 

Y(  113)=Y(  134)*YU:>0)-Y  liiJ)»Y(UV) 

Y(  114)*Y<133)+YU:>o)*YliiO)  +  T'  (  i.50)*Y(  133) 

Y(115)«Y(137)*Y(i:JU)+Y{i^tt)*Y(1^9)+Y(119)*Y(l53)*Y(129)+Y(139)* 
lY(154)+Yi237)*Y(103)*Y(192)-Y(.^3e)*Y(lu3)*Y<191) 
Y(U6)=Y(138)*Y(130)-Y( 137)*Y< 129)+Y< 119)*Y{ 133)*Y( 130)+Y( 140>* 
1Y(155)*Y(238)*Y(105)*Y( 192)+Yi237)*Y< 103>*Y{ 191) 

RETURN 

END 


c 


'FT*' 
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outmuUT Inc  UcRIV( T*C*lO 
TRAJECTORY  PRuokai*i  ouurI 
Cui’irMji't  Y  »'kY 

QIMtNoluN  C  ( 25 )  *u( 25 )  *Y  (oouo )  si£Y  ( 5o )  *2(50) 

Y  ( 108  )  xSORT  ( C ( 4 ) **2+C (  5 ) **2+C { o )  **2 ) 
Y(215)=Y(108)-Y(229) 

CaLL  mkuCFT (Y(215)*Y(l00)*Y(lol)*Y(lo2)*Y(lus ) *Y( 104) ) 
CmLL  Frumka (nY(33)*1*1*Y(215) *uuv  »2 ) 

Y( 102)=Y(102)+Y(102 j#uOV 
Y(  1u'9)=6uRT  (C(5)**2+C(oi**2) 

Y(46)=Y(lu9)/Y(l0o) 

Y(47)*C(4)/Y(  106) 

Y(52)*C(5)/Y( 109) 

Y ( 53 ) *C ( 6 ) / Y ( 109) 

SiiTuP  LL  ma(kIX 
Y( 72)** (46) 

Y ( 73 ) = Y ( 47 ) 

Y  (  7*+  )  =0* 

Y(75)=-Y(h7)*Y(:><1) 

Y(76)*Y(46)*Y(52) 

Y(77)*-Y(53) 

Y(7o)*-YU7)*Y(3j) 

Y(79)=Y(46)*Y(33) 

Y  (  oo  )  *Y  (  5<: ) 

cm;  ScTuP  Or  Ll  ivimTkIX 

CaLL  m\TVtC( Y(72)*C(l)*Y(to©)*0) 

Y(6o)xY(06)-Y(204)*Y< 109) 

I F ( KY ( 7 ) )  9*9*0 

o  CALL  FROriRA  ( \Y  (  34  ) » 1  •  1  •  Y  (  215  )  »  Y  (  227 )  *2  ) 

CALL  FR0nRA<KY(35)*l»l*YUl5)  *Y(22e)  *2) 

CALL  SI  NCOS ( Y ( 22e ) * ow vCw »o  ) 

Y ( 160)*Y(o6)+Y(227)*Cw 
Y ( lol ) aY ( o7 ) 

Y  ( lo2  )  *Y  (  oti  )+Y  ( tt.  7  )  *sw 
00  TO  lu 

9  Y  (  160  )  =  Y  ( t>6  ) 

Y ( 16 1 ) *Y ( 67  ) 

Y( l62)aY(60) 

10  CALL  i*mTVtC<  Y(  72  )  »  Y  ( loo  ) »  Y  <  :)7  )  *1) 

CaLL  i»iaT vc.o(v.(lu)*Y(57)*Y(yo)*4) 

Y(95)=OuRT( Y(90)»*2+Y (9i)~*C+Y (V2)**^) 

Y( ll0)*Y(9S)/( Y ( 1  OS ) * 1 1 16 • 4 ) 

Y< 106 ) =*001168 4S*Y ( 102)*Y(95)**2 
Y(  107)=s2«174*Y(22V)**2/Y(  100 ) **3 
IF( Y ( 2 1 5) -Y (200) ) 15*15*20 
20  Y ( 111 >*0* 

Y ( 1 12 ) =0* 

Y ( 1 1 5 ) B0 • 

Y ( 11h)=0* 

Y ( 1 15 ) *0 • 

Y ( llO)a0# 
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60  TO  16 

15  CALL  aER0( T  tC.U ) 
lo  00  17  1—1*3 

Y  (  I+142)=Y(106)*Y(234)*Y(  1 10  I 

17  Y< 1*145 )=Y( 106 )*Y( 234 )*Y( 235) *Y( I+113) 

CmLL  maTVEC(C( 10 ) *Y ( 143 ) * Y  (  37 ) *0) 

00  io  1=1*5 

v ( I)=Y( 1+36) ✓ Y ( 255) -Y( 107 ) *C  < 1+5) 

10  U( I+5)=C( I ) 

U( 7 ) =(C < 9) *C(85*(Y( 231 )-Y( 252)  )*Y(146)  )/Y(230K- 
A=KY< 12) 

0(8)  =  (C(7)*C(9)*(A*Y(  232)-Y(230)  )4Y(  147)  )/Y(231.* 
U(  9)  =  (C( 7)*C<8 )*( Y( 230)-A*Y(231) )*Y( 148) )/Y<232 1 
0(  10)=U  13)*C(9)-C(  16>*U8) 
0(11)*L(14)*C{9)-C(17)*L(8) 

OC 12)«L( 15)*C(9)-C( lb)*C<8) 

0(13)  M-C  ( 1*0  )  *C  (  9  )  +C  ( 16  )  *U  7  )  *  A 
0( 14)*“C< 11 )*C(9)+C( 17)*C( 7 )*A 
D ( 1 5 ) *-C 1 12 ) *C ( 9 ) *C ( 1 8 ) *C ( 7 ) *A 
D(  16)»C(10)*C<8)-C<  13)*U7)*A 
D(17)«C(11)*C(8)-C(14)*C(7)*A 
0 ( 18 ) =C ( 12 )*C( 8 ) -C ( 15 ) *C ( 7 ) *A 
D ( 19 ) * ( 1 •-A ) *C (7  ) 

00  30  1=1*19 

I F ( AtiS (0(1) ) -10. t-8 >25.25,27 
25  D< I )*0® 

27  IF(ABS(C( !) )-10#E-8)28.28,30 

28  C( I )*0. 

30  CONTINUE 

RETURN 

END 


VOLTS  <4-225 
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SUBROUTINE  TERM( T  »C ,U  *  F ) 

TRAJtCTORY  PROGRAM  6UUR1 
COMMON  Y»K  Y 

DIMENSION  C(25)»D(25),Y(8000)»KY(50) 
I F ( T-Y ( 208 ) ) 7*6*7 

6  C0UNT*0. 

7  F=  1  • 

IF  ( ‘‘V  (  8  )  )  5  » 5 » 1 

1  DO  3  1*1,3 
Y<  300 ) =C  < 10) 

Y ( 301 )=C( 13) 

Y ( 302 ) *C ( 16 ) 

Y ( 303 ) =C ( 1 1 ) 

Y ( 304)*C( 14) 

Y ( 305 ) *C ( 1 7 ) 

Y ( 306 ) *C ( 12  ) 

Y ( 307 ) *C ( 1 5  ) 

Y (  308 ) *C ( 18  ) 

CALL  MATINV(Y(300)»Y(309) *  Y  <  3 1 8  )  ) 

DO  2  J* 10 • 18 

2  C(  J)*(C( J)+Y( J+308)  )  / 2 • 

3  CONTINUE 
COUNT^COUNT+l # 

B  =  K.Y<  10) 

IF(B-COUNT)4,5,5 

4  CALL  MATVEC(C( 10) *C( 10) »Y( 327) .4) 
PRINT  100»<Y< I ) .1*327,335) 

100  FORMAT! IX, 1P9E13.5) 

COUNT*0. 

5  CONTINUE 

I F ( Y ( 2 1 8 ) )  12,12,13 

12  F=Y ( 2 1 5 ) -Y ( 200 ) 

IF  (F)  20,20.13 

13  IF(Y(201 )) 16,16,14 

14  IF ( T-Y ( 201 ) ) 16 , 19 , 1 9 

16  IF(Y(222)-YC202) )50»50»lV 

19  F  =  0, 

20  Y(121)=1.0 
50  RETURN 

END 


VOLTS  64-225 
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SUBROUTINE  OUT ( T ,C*D, ERROR *DMY  »DNYX  *DELT ) 

C  TRAJECTORY  PROGRAM  6UUR1 

COMMON  Y.KY 

DIMENSION  C ( 25  > • D ( 2 5 ) »Y(8000) *KY(50) .ERROR (25) 

10  FORMAT ( lH0*5X»4HTIME*12X*lHH*16X»2HVt*15X»2H\/A»16X*lHM* 

1 14X  * 4HQBAR  *  1 IX  » 7HGAMMA  E/ / ) 

11  FORMAT ( 1H0.4X*4HTIME*10X*6HALPHA-»11X*3HPHI * 10X  *7HPHl  PRI 
1*9X*5HALPHA*10X,4HBETA,13X*1HQ*13X.1HR) 

12  FORMAT ( 1H0*5X*4HTIME,15X,4HDELT* 18X.1HP,19X,4HRTAU*17X* 
14HRPSI *16X*6HPHID0T//) 

13  FORMAT! 1H0*5X*4HTIME,12X,2HFX.15X*2HFy*15X»2HFZ*15X*2HMX* 
115X  *2HMY  .15X.2HMZ//  ) 

14  FORMAT! 1H0 ♦ 5X ,4HT I  ME , 1 1 X *4HDDXR *  1 3X » 4HDDYR *  1 3X ,4HDDZR » 
113X.3HDXR.14X*3HDYR*14X,3HDZR// ) 

1  FORMAT! 3X»4HRUN  * I  3  *95X  *7HF0RMAT  til) 

2  FORMAT! 3Xt I 6 *9 8X , 5HPAGE  *13) 

3  FORMAT! 1PE13.4.1P6E17. 7) 

4  FORMAT! lPE13.4t 1P5E21 . 7 ) 

5  FORMAT 11H0) 

6  FORMAT! 1H1) 

7  FORMAT! 1PE1 3 . 4 t 1P7E 1 5# 5 ) 

IF! T-Y! 208 > ) 105*100*105 

100  KK=0 
N*0 
KZ»0 
J«0 

GO  TO  115 


105 

KK=KK+1 

I F ( Y ( 12 1 )  ) 

110 

*11 

0*115 

110 

IF  (KK-KY! 

9)  ) 

320 

.111*320 

111 

KK  =  0 

I  F ( Y (  142)- 

T)  1 

15* 

320*115 

115 

KZ*KZ+1 

I F ( KZ-46 ) 200 1 1 

5*1 

5 

15 

N*N+1 

DO  30  L  =  1 #  5 

WR I TE ! 6  *  1 ) KY ( 2 ) *L 

WR I TE ( 6  »2 ) KY ( 1  )  tN 

GO  TO  (20t21t22.23t24)  ,L 

20  WR I TE ( 6  1 10 ) 

Ms7*KZ+5000 

WR I TE ( 6 1  3 ) <  Y ! I ) tI*500l#M) 
GO  TO  26 

21  WR I TE ( 6 1 1 1 ) 

M«8*KZ+5315 

WR I TE ( 6 1  7 ) (Y(I ) *  I «5316*M) 
GO  TO  26 

22  WRITE(6tl2) 

M=6*KZ+5675 

WR I TE! 6 14 ) ( Y ( I ) *I*5676*M> 
GO  TO  26 


-V 

4-  T 
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23  WR I TE ( 6 •  13 ) 

M*7*KZ+5945 

WRITEI6.3) (Y(I )  *I*5946.M) 

GO  TO  26 

24  WRITE (6* 14) 

M*7#KZ+6260 

WRIT£(6*3)(Y(I)*I *6261 *M) 

GO  TO  26 
26  WRITE<6*6) 

30  CONTINUE 
J*0 
K*0 
JK*C 
KZ*0 

GO  TO  320 

2 00  Y ( 190 ) *Y ( 204 ) *T 

CALL  SINCO$( Y( 190 ) • SOT  *CQT  •  1 ) 

C  SETUP  LE  MATRIX 
Y ( 37 ) *1 • 

Y <  38 ) *0# 

Y ( 39 ) *0# 

Y ( 40 ) *0# 

Y ( 41 ) *COT 
Y ( 42 ) *-SOT 

Y  <  43 ) *0« 

Y( 44 ) *SOT 

Y  <  45 ) «C0T 

CALL  MATVEC ( Y( 37 ) »C  ( 4 ) *Y(69) #0) 

CALL  ARKTAN*Yi71)*Y<70)*Y<220) *1) 

I F ( ABS( Y ( 220 ) ) “0* 1 )  40*40*41 

40  Y(158)®Y(71)/Y(70) 

Y(94)=Y<229)/Y<70)»Y(71)*< 1*-Y <158 )«*2/3.+Y( 158 )**4/5. ) 
GO  TO  42 

41  Y(94)*Y(229)*Y(220) 

42  Y(190)»SQRT(Y(70)**24Y(71)**2) 

CALL  ARKTAN(Y(69) * Y < 190  )  * Y < 22 1 )  *  1 ) 

IF ( A0S( Y ( 221 ) ) -Oel )  43*43,44 

43  Y(159)*Y(69)/Y(190) 

Y<95)«Y(229)/Y<190)*Y(69)*<l*-Y<1595**2/3.+Y<159>**4/5.) 
GO  TO  45 

44  Y(95)*Y(229)*Y(221) 

'.5  Y  ( 141 )  «  (  «  Y  C  92  J  *  C  Y «  144 ) /Y < 110 ) -C ( 9 ) *Y < 90 ))-Y< 91 )* 
1(Y(145)/Y(110)+C(8)*Y(90> ) ) / ( Y < 92 ) **2+Y ( 91 ) **2 ) )+C< 7 ) 

Y  ( 190 ) *SQRT ( Y ( 66 ) **2+Y ( 68 ) **2 ) 

CALL  ARKTAN( Y( 67 ) ,Y(190)*Y(218),0) 

Y ( 2 16 ) *SQRT ( Y ( 66 ) **2+Y  <  67 ) **2+Y ( 68 ) **2  ) 

IF  (J)  310*300*310 
300  J  *  1 
K*l 
JK«1 

GO  TO  311 
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TABLK  6  (CoatlBMd) 


310  J* J+7 
K*K+8 
JK*JK+6 

311  Y ( J+5000 ) =T 
Y(J+5001)SY(215) 
Y< J+5002)=Y(216) 
Y< J+5003)»Y<93) 

Y< J+5004)*Y( 110) 
Y(J+5005)SY(106) 
Y(J+5006)*Y(218) 
Y(K+5315)=T 

Y ( K+5316 ) *Y( 222 ) 

Y ( K+5317 ) =Y ( 223 ) 

Y ( K+5318 ) *Y ( 151 ) 

Y ( K+5319 ) =Y ( 299 ) 

Y ( K+5320 ) *Y ( 298 ) 

Y ( K+532 1 ) =C ( 8 ) 
Y(!C+5322)=C(9) 

Y ( JK+5675 ) *T 
Y( JK+5676)*DELT 
Y(JK+5677)*C(7 ) 

Y( JK+5678)*Y<94) 
Y( JK+5679)*Y<95) 
Y< JK+5680)«Y(141) 


Y( J+5945 

)  *  T 

Y( J+5946 

)  =Y  ( 143 ) 

Y ( J+5947 

)  *  Y  ( 144 ) 

Y( J+5948 

)*Y( 145) 

Y ( J+5949 

)*Y( 146) 

Y( J+S950 

)*Y( 147) 

Y ( J+595 1 

) * Y ( 148  ) 

Y  ( J+6260 

)  *  T 

Y( J+6261 

)  =0  (  1 ) 

Y ( J+6262 

)  *D  (  2  ) 

Y( J+6263 

)  ED  (  3  ) 

Y(  J+6264 

)  =D  {  4  ) 

Y ( J+6265 

) =D { 5  ) 

Y ( J+6266 

)=U( 6) 

IF( Y( 121 

) ) 315  *315*  15 

315 

I F ( KZ-45 

)  320*15,320 

320 

Y ( 142 )*T 

350 

RETURN 

END 


worn  14-MI 
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SUBROUTINE  FN0L3( J*N*G*L*M* ENE.X.Y.D, DERI V.TERM, OUTPUT) 

C  FORTRAN  IV  11-63 

DIMENSION  Y{50)*D(50) *YB( 30*6) »GI 2 (30) >613(30) tGI 4(30) #EFJ30)* 
1EFK30) *EF2(30i *EF3(30) *Yl(30) .ERROR ( 30) *HA(30) *YA( 50 I *OA(5o> * 
2YC(30)*YP(30) *YD(50) 

DOUBLE  PRECISION  XD*YD*YA*YC*YP.Y1 
EC*Y  ( N+3 ) 

1  H*6 

2  HZ*H 

3  LN*N+MAX0(  L*3 ) 

4  NA«0 

5  NB*1 

6  NF*0 

7  NG*0 

8  F*0* 

9  FA*0* 

10  FB=0. 

11  FC*0* 

12  FD*0. 

13  NE*ENE 

DO  200  1*1. LN 
200  YD ( I ) *DBLE ( Y ( I  )  ) 

XD*DBLE(X) 

14  IF ( J-3 )15»21»15 

15  IF(NE) 18*16*18 

16  JA*4 

17  GO  TO  22 

18  RE1»10.**(-ENE) 

19  RE2*l0.**(-ENE-3.0) 

20  REM*10* ** ( -ENE-1 *5 ) 

21  JA*  1 

22  DO  25  I  *  1  * N 

23  DO  24  I C* 1  * 5 

24  YB(  1*10*0. 

25  ERROR ( I ) *0. 

26  CALL  DERIV(X.Y*D) 

DO  300  1*1. N 

G 1 2  < I )*D( I ) 

G13(I )*D(I  > 

GI 4 ( I )*D< I ) 

300  EF( I )*D(  I ) 

27  CALL  OUTPUT (X*Y*D*ERROR*N*L*H) 

28  FD* Y ( N+l ) 

29  I F ( J-2 )  30*129*30 

30  GO  70(31*37*35. 37 )*JA 

31  DO  33  I *1 *LN 

32  YA ( I ) *YD( I ) 

33  DA ( I ) *D (  I  ) 

34  GO  TO  37 

35  HB*H 

36  H*2.*H 
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37  HD2  =  • 5*H 

DO  39  I *1 *N 

38  YB< I ,NB )  *D  {  I ) 

XL  =  D( I )  *  HD 2 

39  Y { I ) *SNGL ( YD ( I )+XL) 

X=SNGL ( XD+HD2 ) 

40  CALL  DERI V  (X.Y.GI2 ) 

41  DO  42  I  *1  »N 

XL  »  G I  2 ( I ) *HD2 

42  Y ( I ) =  SNGL ( YD( I )+XL) 

43  CALL  DERI V  (X,Y*GI3) 

44  DO  45  I =1 » N 
XL*G I  3 ( I )*H 

45  Y< I )=SNGL( YD( I )+XL) 

X*SNGL ( XD+H ) 

46  CALL  DERIV(X»Y»GI4) 

47  HD6  *H/6* 

GO  TO(48.55f60#66)»JA 

48  DO  52  I =1 *N 

XL* ( D ( I )  +  2.*(GI2(I)  +  G I  3 ( I ) )  +GI4(I)J*HD6 

49  YC( I )=YD( I J+XL 

51  YD { I )=YA(  I  ) 

52  ERROR ( I ) =0 • 

53  JA*3 

54  GO  TO  35 

55  DO  57  1=1, N 

XL* ( D ( I  )  +  2**(GI2(  I  )  +  GI 3 (  I  )  )  +GI4(I))*HD6 

56  YD { I ) =YD ( I ) +XL 

57  ERROR ( I )*SNGU YD( I) -YP ( I ) ) / 15. 

58  JA= 1 

59  GO  TO  681 

60  DO  62  1=1, N 

61  YD  C I  I =YC ( I  ) 

XL  * ( D (  I  )  +  2«* ( GI 2 (  I  )  +  G I  3 (  I  )  )  +GI4(I))*HD6 

62  YP ( I ) =YA { I J+XL 

63  H*HB 

64  JA*2 

65  GO  TO  681 

66  DO  68  I  =1 , N 

XL* ( D ( I )  +  2a*(GI2(I)  ♦  G I  3 ( I )  )  +GI4(I))*HD6 

67  YD< I )=YD( I ) +XL 

68  ERROR ( I ) *0 a 
681  DO  69  1=1,  N 

69  Y (  I ) =SNGL ( YD( I ) ) 

XD*XD+H 

X*SNGL I XD ) 

70  CALL  DERI V CX,Y,D) 

71  FC*F 

72  CALL  TERM(X»Y, D,F) 

73  I F ( ABS ( F ) -1 aOE-5  >731,731,733 
731  NF*5 
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732  GO  TO  124 

733  IF  <  F ) 74, 124,76 

74  FA*1 • 

75  GO  TO  77 

76  FB*1. 

77  IF(FA-FB>83, 78*83 

78  NF=NF+1 

79  JA*4 

80  NB*1 

81  H*H*F/(FC-F) 

82  I F ( NF-4 '37,37,124 

83  I F C NE)84*117*84 

84  I F ( JA-1 ) 117*85*117 

85  IF(J-3)86,117*86 

86  DO  95  1*1, N 

I F ( Y  C I ) )886,885,886 

885  HA  (  I ) *1000, 

GO  TO  95 

886  IF( EC) 880*890*87 

87  I F( ABS< Y( I ) )-EC )  880*880*890 

880  IF(  ABS(ERRORU)  >-RE2>  882,94,881 

881  1F<ABS<ERROR<  I ) >-REl > 94,94,882 

882  HA( n*H*< REM/ <ABS<ERROR<  I ) > +.0000000001 ) >** <  *2 ) 

883  GO  TO  95 

890  IF<ABS<ERROR<  I )/Y< I ) )-RE2 ) 892.94,891 

891  I F  < AbS< ERROR < I )/Y< I ) )-REl ) 94, 94, 892 

892  HA( I )*H*< REM/ <ABS<ERROR<  I ) /Y< I) )♦• 0000000001 >  >**< #2) 

893  GO  TO  95 

94  HA ( I ) *H 

95  CONTINUE 

96  HB«HA(N) 

97  DO  98  1*1, N 

98  HB*AMIN1 <HA< I ) ,HB) 

99  IF (H-HB) 100*117*101 

100  IF < HZ-H ) 101,101*116 

101  DO  103  1*1, LN 

102  YD ( I ) *YA ( I ) 

Y< I ) *SNGL ( YD  < I ) ) 

103  D< I ) *DA ( I ) 

104  IF(NB-6)107, 105,105 

105  XD*XD-H 

106  GO  TO  109 

107  XD=XD-2.*H 

108  HZ*H 

109  H*HB 
X*SNGL(XD) 

CALL  DERI V(X,Y »D) 

110  NB= 1 

111  XABS*ABS( •000001*X) 

112  I F ( ABS ( H ) -XABS ) 113,113*117 

113  NG*NG+1 
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114  H*SIGN( XABS.HB) 

115  IFiNG-10) 124*126.126 

116  HZ*H  , 

117  IFCM) 118*118*121 

118  IF(  ABS{  Y(N-*-l  )-FD)-Y(N+2  )) 29, 119, 119 

119  FD*Y(N+1) 

120  GO  TO  124 

121  NA=NA+1 

122  I F ( M-NA *123,123,29 

123  NA*0 

124  CALL  OUTPUT(X*Y»D,ERROR,N,L*H) 

125  lF(NF-4)29, 29,126 

126  WRITE  (6.127) 

127  FORMAT ( 1H0 ) 

128  RETURN 

129  NB*NB+1 

130  I F ( NB-6 )30#131*136 

131  DO  134  1  =  1, N 

132  EF3 ( I )  *YB ( 1,3) 

133  EF2  < I ) * YB ( I *4  > 

134  EF1 ( I ) *YB ( 1,5) 

135  GO  TO  137 

136  NB= 10 

137  HD24  =H/24* 

DO  138  1  =  1, N 

XL  =(55.*D(I)  -59,*EF1 ( I )  ♦37.*EF2(I)  -9,*EF3 ( I ) ) *HD24 
YP( I )*YD< I )4XL 

138  Y ( I ) *SNGL ( YP ( I  )  ) 

X*SNGL( XD+H) 

139  CALL  DERIV(X,Y,EF) 

140  DO  142  1=1, LN 

141  YA( I )=YD( I  ) 

142  DA ( I ) =D (  I  ) 

143  DO  148  1  =  1, N 

XL  *  (9.#EF(I)  +19,*D( I )  -5,*EF1(I)  ^EF2 ( I ) )*HD24 

144  YD  (  I  )  *YD  (  I  )  -f  XL 

145  ERROR ( I ) =-5NGL ( YD { I ) -YP < I ) )/14, 

146  EF3 ( I ) »EF2 ( I ) 

147  EF2 ( I )*EF1 ( I ) 

148  EF1 ( I ) >D( I ) 

149  GO  TO  681 
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TABLE  8 


KY  ARRAY  VARIABLES 


KY  (1) 

Date 

rar(3S)  aw 

KT(3<!)-Klf(50) 

KY  (2) 
KY  (3) 

Run  No. 

J 

KY  (4) 

N 

KY  (5) 

L 

KY  (6) 

M 

KY(7) 

Wind 

KY  (8) 

L^Bj  correction 

KY  (9) 

Print  Freq. 

WT 

KY (10) 

Print  Freq.  j  l  "] 

L  Bj 

KY  (11) 

Aero  Option 

KY(12) 

k 

KY(13)-KY(20)  Available  for  input 
KY (21)  C 

x 

KY (22)  C 

x 

P 

KY (23)  C 

y 

KY (24)  C 

z 

KY (25) 

KY (26)  C 

6 

KY (27)  Ct 

P 

KY  (28)  C 

m 

KY  (29)  C 

n 

KY  (30)  C 

in 

q 

KY  (31)  C 

n 

P 

KY (32)  C 

TO. 

6a 

KY  (33)  D 
KY (34)  Vw 


Available  for 
input 


» 


i 


,5 

t  i 
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TABLE  9 


Y  ARRAY 

VARIABLES 

Y(l)  - 

Y(36) 

Setup  only 

Y  (37) 

-  Y(65) 

Intermediate  storage 

Y  (66) 

V 

XLE 

Y  (67) 

V 

yLE 

Y  (68) 

V 

ZLE 

Y  (69) 

*E 

Y  (70) 

YE 

Y  (71) 

ZE 

Y  (72) 

-  Y  (80)  j 

'lh J  array 

Y  (81) 

-  Y  (89) 

Available  for 

program 

Y  (90) 

V 

XBA 

Y  (91) 

V 

yBA 

Y  (92) 

V 

ZBA. 

Y  (93) 

VA 

Y  (94) 

R 

T 

Y  (95) 

R. 

0 

Y  (96) 

-  Y(104) 

Intermediate 

storage 

Y(105) 

C 
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Y  ARRAY  VARIABLES  (Continued) 

Y (106)  Q 

Y (107)  -  Y(109)  Intermediate  Storage 

Y (110)  U 

Y(lll)  -  Y(116)  Intermediate  Storage 

Y (117)  Cx 

P 

Y(I18)  C, 

Y (119)  Cn 

P 

Y (120)  -  Y (130)  Intermediate  Storage 

Y(131)  0' 

Y (132)  Cx 

Y (133)  C 

Y (134)  C 

z 

Y (135)  v 

Y (136) 

P 

Y(I37)  C 

m 

Y(138)  Cn 

Y(139)  C  (pitch) 

fflq 

Y(I4o)  C  (ya v) 

m 

q 

Y (141)  b 

Y(142)  Intermediate  Storage 
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Y  ARRAY  VARIABLES  (Continued) 

Y (143)  Fx 

Y(144)  Fy 

Y(145)  F„ 

z 

Y(146)  Mx 

Y (147)  M 

y 

Y (148)  M 

z 

Y (149)  -  Y (162)  Intermediate  Storage 

Y(163)  -  Y(189)  Setup  only 

Y(190)  -  Y (199)  Intermediate  Storage 

Y (200)  ht  (ft) 

Y (201)  tt  (sec) 

Y (202)  at  (deg) 

Y (203)  Ag  (ft) 

Y (204)  (rad/sec) 

Y  (205)  Error 

Y(20&)  ENE 

Y (207)  At  (sec) 
o 

Y (208)  t  (sec) 

o 

Y(209)  -  Y (214)  Available  for  input  data 

Y  (215)  h  (ft) 

Y (216)  VE  (ft/sec) 


NOLTR  64-225 


Y  ARRAY  VARIABLES  (Continued) 


Y  (217) 

ae 

(deg) 

Y (218) 

7E 

(deg) 

Y  (219) 

9 

(deg) 

Y (220) 

TE 

(deg) 

Y  (221) 

(deg) 

Y (222) 

5 

(deg) 

Y  (223) 

0 

(deg) 

Y (224) 

p 

(rad/sec) 

Y  (225) 

q 

(rad/sec) 

Y  (226) 

r 

(rad/sec) 

Y (227) 

vw 

(ft/sec) 

Y  (228) 

Aw 

(deg) 

Y (229) 

R_ 

E 

(ft) 

Y (230) 

I 

X 

(slug— ft2 ) 

Y  (231) 

I 

y 

(slug-ft3) 

Y  (232) 

i 

z 

(slug-ft*) 

Y  (233) 

in 

(slug) 

Y  (234) 

A 

(ft*) 

Y (235) 

d 

(ft) 

Y  (236) 

6 

(rad) 

Y (237) 

5a 

(rad  or  deg  depending  on  C 
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Y  ARRAY  VARIABLES  (Continued) 


Y (238)  6_  (rad  or  deg  depending  on  C  ) 

D  IB. 

*A 

Y(239)  -  Y(297)  Available  for  input 
Y (298)  0 

Y (299)  a 

Y (300)  -  Y(335)  Intermediate  storage 
Y(336)  -  Y(599)  Available  for  program 
Y (600)  -  Y (4999)  Table  storage 
Y(5000)  -  Y(8000)  Output  storage 
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TABLE  10 

c 

ARRAY  VARIABLES 

C(l) 

*R 

C(12) 

Si 

C(2) 

C(13) 

^12 

C(3) 

*R 

C(14) 

l22 

C(4) 

XR 

C(15) 

l32 

C(5) 

yr 

C(16) 

ll3 

C(6) 

ZR 

C(17) 

l23 

C(7) 

P 

C(18) 

l33 

C(8) 

q 

C(19) 

<f>' 

C(9) 

C(X0) 

r 

ln 

C(20) 

C(21) 

Blank 

C(ll) 

''21 

C(22) 

Error  Control 
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TABLE  11 

D  ARRAY  VARIABLES 


D(l)  XR 

D(2)  ¥r 

D(3)  2r 

D(4)  XR 

D(5)  tR 

D(6)  ZR 

D<7)  p 

D(8)  q 

D(9)  r 


D(10)  lu 

D(ll)  l21 

D(12)  l31 

D(13)  112 

0(14)  J'22 

0(15)  i32 

D(16)  t13 

0(17)  S3 

0(18)  *33 

D(19) 
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FORMAT 

t 

(sec) 

FORMAT 

t 

(sec) 

FORMAT 

t 

(sec) 

FORMAT 

t 

(sec) 

FORMAT 

t 

(sec) 


TABLE  12 


OUTPUT  FORMAT 


h  E  A 

(ft)  (ft/sec)  (ft/sec) 


M  Q  'E 

(lb/ft*)  (deg) 


5  0  Q  a  p  q  r 

(deg)  (deg)  (deg)  (deg)  (deg)  (rad/sec)  (rad/sec) 


At 

(sec) 


(rad/sec) 


Rr  r*  ft 
(ft)  (ft)  (rad/sec) 


F  F  F  M  M  U 

x  y  z  x  y  z 

(lb)  (lb)  (lb)  (ft-lb)  (ft-lb)  (ft-lb) 


^R  ^R  2R  *R  *R  *B 

(ft/sec*)  (ft/sec*)  (ft/sec*)  (ft/sec)  (ft/sec)  (ft/sec) 
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F,G-  I  AERODYNAMIC  DATA 


and  BODY  AXES  SYSTEMS 
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FIG. 2  INERTIAL  AND  LOCAL  AXES  SYSTEMS 
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*r,Xe 


FIG. 5  EARTH  AND  INERTIA1.  AXES  SYSTEMS 
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